#ifndef MLTool_TRIANGULAR_MATRIX_H
#define MLTool_TRIANGULAR_MATRIX_H

#include"layout.h"
#include"vector.h"
#include"tags.h"
#include<functional>
#include<utility>
#include<cassert>

namespace math{
template<class T, class TRI>
class triangular_matrix : public vector<T>
{
public:
    typedef vector<T> Array;
    typedef typename Array::difference_type difference_type;
    typedef typename Array::size_type size_type;
    typedef typename Array::value_type value_type;
    typedef typename Array::pointer pointer;
    typedef typename Array::const_pointer const_pointer;
    typedef typename Array::reference reference;
    typedef typename Array::const_reference const_reference;
    typedef typename Array::iterator iterator;
    typedef typename Array::const_iterator const_iterator;
    typedef typename Array::reverse_iterator reverse_iterator;
    typedef typename Array::const_reverse_iterator const_reverse_iterator;
    typedef matrix_tag category;
    typedef typename TRI::triangular_category triangular_category;

    using vector<T>::data_;


    template<class U, class TRI2>
    friend class triangular_matrix;

    triangular_matrix(size_type i = 0, size_type j = 0):
        rows(i), columns(j), vector<T>(TRI::storage_size(i, j))
    {}


    triangular_matrix(size_type i, size_type j, const T &init):
        rows(i), columns(j), vector<T>(TRI::storage_size(i, j), init)
    {}


    template<class U>
    triangular_matrix(const triangular_matrix<U, TRI> &m)
        :rows(m.rows), columns(m.columns),
          vector<T>(m)
    {}


    //swap
    void swap(triangular_matrix & m)
    {
        std::swap(rows, m.rows);
        std::swap(columns, m.columns);
        vector<T>::swap(m);
    }

    inline friend void swap(triangular_matrix & m1, triangular_matrix & m2)
    {
        m1.swap(m2);
    }

    inline size_type size1() const{return rows;}
    inline size_type size2() const{return columns;}
   
    inline T& operator()(size_type i, size_type j)
    {
        return data_[TRI::element(i,j,rows,columns)];
    }

    inline const T& operator()(size_type i, size_type j) const
    {
        if(TRI::zero(i,j))
            return zero_;
        else if(TRI::one(i,j))
            return one_;
        else
            return data_[TRI::element(i,j,rows,columns)];
    }

	//operator +
    friend inline triangular_matrix operator +(const triangular_matrix &m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "1st called!" << std::endl;
        triangular_matrix ans(m1.rows, m1.columns);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] + m2.data_[i];
        return ans;
    }

    friend inline triangular_matrix&& operator +(triangular_matrix &&m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] += m2.data_[i];
        return std::move(m1);
    }

    friend inline triangular_matrix&& operator +(const triangular_matrix &m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] += m1.data_[i];
        return std::move(m2);
    }

    friend inline triangular_matrix&& operator +(triangular_matrix &&m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] += m2.data_[i];
        return std::move(m1);
    }

	//operator -
    friend inline triangular_matrix operator -(const triangular_matrix &m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "1st called!" << std::endl;
        triangular_matrix ans(m1.rows, m1.columns);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] - m2.data_[i];
        return ans;
    }

    friend inline triangular_matrix&& operator -(triangular_matrix &&m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] -= m2.data_[i];
        return std::move(m1);
    }

    friend inline triangular_matrix&& operator -(const triangular_matrix &m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] = m1.data_[i] - m2.data_[i];
        return std::move(m2);
    }

    friend inline triangular_matrix&& operator -(triangular_matrix &&m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] -= m2.data_[i];
        return std::move(m1);
    }

	//operator *
    friend inline triangular_matrix operator *(const triangular_matrix &m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "1st called!" << std::endl;
        triangular_matrix ans(m1.rows, m1.columns);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] * m2.data_[i];
        return ans;
    }

    friend inline triangular_matrix&& operator *(triangular_matrix &&m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] *= m2.data_[i];
        return std::move(m1);
    }

    friend inline triangular_matrix&& operator *(const triangular_matrix &m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] *= m1.data_[i];
        return std::move(m2);
    }

    friend inline triangular_matrix&& operator *(triangular_matrix &&m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] *= m2.data_[i];
        return std::move(m1);
    }

	//operator /
    friend inline triangular_matrix operator /(const triangular_matrix &m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "1st called!" << std::endl;
        triangular_matrix ans(m1.rows, m1.columns);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] / m2.data_[i];
        return ans;
    }

    friend inline triangular_matrix&& operator /(triangular_matrix &&m1, const triangular_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] /= m2.data_[i];
        return std::move(m1);
    }

    friend inline triangular_matrix&& operator /(const triangular_matrix &m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] = m1.data_[i] / m2.data_[i];
        return std::move(m2);
    }

    friend inline triangular_matrix&& operator /(triangular_matrix &&m1, triangular_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.rows == m2.rows);
       assert(m1.columns == m2.columns);
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] /= m2.data_[i];
        return std::move(m1);
    }

 	//operator - (负号)
	friend inline triangular_matrix operator -(const triangular_matrix &m)
	{
		triangular_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = -m.data_[i];
		return ans;
	}

	friend inline triangular_matrix&& operator -(triangular_matrix &&m)
	{
		m *= -1;
		return std::move(m);
	}
    // 矩阵和标量的运算
    //operator +
	friend inline triangular_matrix operator +(const triangular_matrix &m, value_type t)
	{
		triangular_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] + t;
		return ans;
	}

	friend inline triangular_matrix operator +(value_type t, const triangular_matrix &m)
	{
		return m + t;
	}

	friend inline triangular_matrix&& operator +(triangular_matrix &&m, value_type t)
	{
		m += t;
		return std::move(m);
	}


	friend inline triangular_matrix&& operator +(value_type t, triangular_matrix&&m)

	{
		m += t;
		return std::move(m);
	}

    // operator -
	friend inline triangular_matrix operator -(const triangular_matrix &m, value_type t)
	{
		triangular_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] - t;
		return ans;
	}

	friend inline triangular_matrix operator -(value_type t, const triangular_matrix &m)

	{
		triangular_matrix ans(m.size1(). m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = t - m.data_[i];
		return ans;
	}


	friend inline triangular_matrix&& operator -(triangular_matrix &&m, value_type t)

	{
		m -= t;
		return std::move(m);
	}


	friend inline triangular_matrix&& operator -(value_type t, triangular_matrix&&m)

	{
		for(size_type i = 0; i < m.size(); ++i)
			m.data_[i] = t - m.data_[i];
		return std::move(m);
	}
	
	// operator *
	friend inline triangular_matrix operator *(const triangular_matrix &m, value_type t)
	{
		triangular_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] * t;
		return ans;
	}

	friend inline triangular_matrix operator *(value_type t, const triangular_matrix &m)
	{
		return m * t;
	}

	friend inline triangular_matrix&& operator *(triangular_matrix &&m, value_type t)
	{
		m *= t;
		return std::move(m);
	}

	friend inline triangular_matrix&& operator *(value_type t, triangular_matrix &&m)

	{
		m *= t;
		return std::move(m);
	}

	//operator /
	friend inline triangular_matrix operator /(const triangular_matrix &m, value_type t)
	{
		triangular_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] / t;
		return ans;
	}

	friend inline triangular_matrix operator /(value_type t, const triangular_matrix &m)
	{
		triangular_matrix ans(m.size1(). m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = t / m.data_[i];
		return ans;
	}

	friend inline triangular_matrix&& operator /(triangular_matrix &&m, value_type t)

	{
		m /= t ;
		return std::move(m);
	}

	friend inline triangular_matrix&& operator /(value_type t, triangular_matrix &&m)

	{
		for(size_type i = 0; i < m.size(); ++i)
			m.data_[i] = t / m.data_[i];
		return std::move(m);
	}

private:
    size_type rows;
    size_type columns;
    static const T zero_;
    static const T one_;
};

template<class T, class TRI>
const T triangular_matrix<T, TRI>::zero_ = T();

template<class T, class TRI>
const T triangular_matrix<T, TRI>::one_ = static_cast<T>(1);

}
#endif
